% ESTIMA_INITIAL 
% loading starting parameter values for estimation routine
% called by ESTIMA
%
%
% For: Alejandro Justiniano's LMJ Project (Spring 2010)
% RA: Xi Luo
% File Created on 04/08/2010
prpar.prior = prpar_mat;
prpar.alphap = prpar_alpha;
prpar.betap = prpar_beta;
prpar.lbnd = prpar_lb;
prpar.ubnd = prpar_ub;

if flag_prssload==1
prss.prior = prss_mat;
prss.alphap = prss_alpha;
prss.betap = prss_beta;
prss.lbnd = prss_lb;
prss.ubnd = prss_ub;
else 
    prss=[]; 
end 
% 
% tryme=genmultrand(tempdist,prpar_alpha(parposest),prpar_beta(parposest),1000);
% nsol=0;
% for ii=1:size(tryme,1);
%     temp=parloadval;
%     temp(parposest)=tryme(ii,:)';
%     [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag]=modoriglin(temp,solveopt,addsol);
%     if isequal(eu,[1;1]);
%         nsol=nsol+1;
%         lpostd =lmj_postmax( temp(parposest),parloadval, ...
%                 parposest,funcmod,Y,trainvec,prpar,prss, solveopt, addsol, ssposest); 
%     end
% end
%             
%%            
display('________________________________'); 
display('Begin Generating Starting Values'); 

npoints_full = 0; 
if case_initial == 1
    if flag_strvalsload == 0 
        error('case_intial = 1; but starting parameter values not loaded')
    else 
        disp('case_initial = 1, and starting parameter values have been loaded')
        npoints_full = 1; 
    end
    nstrvals=size(parstrvals,2); 
elseif case_initial == 2 || case_initial == 3
    parstrvals = zeros(numpar, nstrvals); 
    icount = 0; % iteration counter
    fcount = 0; % fill count; 
    offset=1e-7;
    rand('twister',seednum);
    disp('_________________________________________________');
    disp('Starting values, uniform from grid');
    pause(0.25);
    
    disp('_________________________________________________');
    disp('Starting values, uniform from prior');
    pause(0.25);
    
    while npoints_full == 0; 
        admissible = 0; 
        icount = icount + 1; 
        parstrval_temp = zeros(numpar, 1);
        switch case_initial 
            case 2 % random from uniform 
              
                for ii = 1:numpar; 
                    if  pargrid(ii, 2) - pargrid(ii, 1) > 2*offset 
                        parstrval_temp(ii) = unifrnd(pargrid(ii, 1)+offset, pargrid(ii, 2)-offset); 
                    else 
                        parstrval_temp(ii) = unifrnd(pargrid(ii, 1), pargrid(ii, 2));
                    end     
                end
                parstrval_temp(parposcal) = parcalval; 
        case 3 
            % random from prior 
            tempdist = char(prpar_dist(parposest));                                
            parstrval_est_temp = genmultrand(tempdist,prpar_alpha(parposest),prpar_beta(parposest),1);
            parstrval_temp(parposest) = parstrval_est_temp; 
            parstrval_temp(parposcal) = parcalval; 
        end
        
        tone=parstrval_temp(parposest)-prpar_lb(parposest);
        two=prpar_ub(parposest) - parstrval_temp(parposest); 
        if all(tone >= 0) == 1 && all(two >= 0) == 1; 
            % Evaluate likelihood and posterior 
           [lpostd,likel,lpriord,lprior_ss,flag_ok]=feval(@lmj_postmax,parstrval_temp(parposest),parstrval_temp,...
                parposest, funcmod, Y, trainvec, prpar, prss, solveopt, addsol, ssposest); 
            
            disp(' '); 
            disp(sprintf('LogLikel=%10.2f',likel));
            disp(sprintf('lpriord=%10.2f',lpriord));
            disp(sprintf('lprior_ss=%10.2f',lprior_ss));        
            disp(sprintf('LogPosterior=%5.5e',lpostd)); 
            
            if lpostd > -1e15 && flag_ok==1 
                admissible = 1; 
            end 
        end 
            if admissible == 1;
                fcount = fcount + 1; % Keep or not
                parstrvals(:, fcount) = parstrval_temp; 
            end
            npoints_full = (fcount == nstrvals);  
    end 
else 
    error('varialbe CASE_INITIAL only takes integer values on 1, 2 and 3');
end
%%
disp('Checking initial densities....'); 
logprior_mat = zeros(nstrvals, 2);
logpostd_mat  = zeros(nstrvals, 2); 
loglikel_mat = zeros(nstrvals, 2); 
for ii = 1:nstrvals;
    parvec = parstrvals(:, ii); 
    parest = parstrvals(parposest, ii); 
    
    logpostd_mat(ii, 1) = ...
        lmj_postmax(parest, parvec, parposest, funcmod, Y, trainvec, prpar, prss, solveopt, addsol, ssposest); 
    
    if logpostd_mat(ii,1) < -1e15 
        error('Likelihood too small! Parameter value is not admissible') 
    end 
            
    loglikel_mat(ii, 1) = ...
        lmj_likelcentral(parest, parvec, parposest, funcmod, Y, trainvec, solveopt, addsol);    
    dispaj('Starting densities for candidate ',ii); 
    disp([sprintf('LogPosterior=%5.5e',logpostd_mat(ii, 1)),sprintf(' LogLikelihood=%5.5e',loglikel_mat(ii, 1))]); 

    % sum of prior densities on parameters and on steady state variables
    logprior_mat(ii, 1) = logpostd_mat(ii, 1) - loglikel_mat(ii, 1);
    
end 
disp('Obtained PARSTRVALS'); 
savecell([parnames num2cprec(parstrvals,7)],[],outpath,'starting values'); 

